For this analysis I use the HCC data set with the relative water levels (counterclimatic and observed). Data can be accessed here
In this part we will make a plot of the regional observed water level changes and compare it to the global average. What we use for this is the hourly hcc waterlevel data that can be downloaded for every year seperately. In order to get a first idea of this data, we can take one example year: For this we extract the relevant variables and then take Stockholm as location.
setwd("/Users/hannascheiwe/Desktop/CLEWS2.0/uni_kram/CC_Adaptation/data/hourly_data_relative_waterlevels_observed")
wl1997 <- nc_open("hcc_obsclim_waterlevel_global_hourly_1997.nc")
print(wl1997)## File hcc_obsclim_waterlevel_global_hourly_1997.nc (NC_FORMAT_NETCDF4):
##
## 3 variables (excluding dimension variables):
## double waterlevel[stations,time] (Contiguous storage)
## comment: This variable represents the relative coastal water level as a combination of reconstructed geocentric water level data based on Dangendorf et al. (2019), vertical land motion (VLM) based on Oelsmann et al. (2023) and hourly storm-tide variations based on Muis et al. (2020). The data is aligned with the geocentric water level such that the last time step in the record is identical. It therefore describes the water level relative to a virtual tide gauge that is at the height of the geocentric water level at the last date in the record.
## units: mm
## positive: up
## coordinates: lat lon
## long_name: Relative Water Level
## double lat[stations] (Contiguous storage)
## standard_name: latitude
## long_name: Latitude
## units: degrees_north
## axis: Y
## double lon[stations] (Contiguous storage)
## standard_name: longitude
## long_name: Longitude
## units: degrees_east
## axis: X
##
## 2 dimensions:
## stations Size:18719 (no dimvar)
## time Size:8760
## long_name: Time
## axis: T
## units: hours since 1979-01-01
## calendar: proleptic_gregorian
##
## 9 global attributes:
## Convention: CF-1.9
## featureType: timeSeries
## references: The dataset is described in Treu et al. (2023) <https://doi.org/10.5194/essd-2023-112>. The source datasets are described in Dangendorf et al. (2019) <https://doi.org/10.1038/s41558-019-0531-8>, Muis et al. (2020) <https://doi.org/10.3389/fmars.2020.00263> and Oelsmann et al. (2023) <https://doi.org/10.5281/zenodo.8308347>
## institution: Potsdam Institute for Climate Impact Research (PIK)
## project: Inter-Sectoral Impact Model Intercomparison Project phase 3a (ISIMIP3a)
## contact: ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>
## isimip_id: d193c353-beee-4e8e-bdd4-84d7a6b74297
## grid_info: The irregular coastal grid points are distributed equidistantly along a smoothed global coastlines with a distance of 50 km globally and 10km in Europe with additional points added for tide gauge stations (described in Muis et all. 2019, https://doi.org/10.3389/fmars.2020.00263).
## title: Hourly coastal relative water levels
water_level_1997 <- ncvar_get(wl1997, "waterlevel")
lat <- ncvar_get(wl1997, "lat")
lon <- ncvar_get(wl1997, "lon")
time_values_1997 <- ncvar_get(wl1997, "time")
stockholm_lat <- 59.3293
stockholm_lon <- 18.0686
distances_1997 <- sqrt((lat - stockholm_lat)^2 + (lon - stockholm_lon)^2)
closest_station_idx_1997 <- which.min(distances_1997)
stockholm_water_level_1997 <- water_level_1997[closest_station_idx_1997 , ]
# Calculate global average water level for each time step (average across all stations)
global_average_water_level_1997 <- rowMeans(water_level_1997, na.rm = TRUE)
mean_water_level_all_stations_1997 <- apply(water_level_1997, 2, mean, na.rm = TRUE)
#--> calculates the average water level across all stations for each time step
global_avg_sea_level_1997 <- mean(mean_water_level_all_stations_1997, na.rm = TRUE)Okay now we do have the average values we were looking for, but the goal is to have a time series over the whole period. So the next step will be to get the mean value of all stations for the respective year in order to get the development of the global average sea level change. For this we can use a loop that applies all these steps for the different data. Here it is important that all the files are stored in the same folder.
global_avg_sea_levels <- vector("numeric", length = 19)
stockholm_water_levels <- vector("list", length = 19)
for (year in 1997:2015) {
file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
if (file.exists(file_name)) {
dat <- nc_open(file_name)
water_level <- ncvar_get(dat, "waterlevel")
lat <- ncvar_get(dat, "lat")
lon <- ncvar_get(dat, "lon")
time_values <- ncvar_get(dat, "time")
distances <- sqrt((lat - stockholm_lat)^2 + (lon - stockholm_lon)^2)
closest_station_idx <- which.min(distances)
stockholm_water_level <- water_level[closest_station_idx , ]
stockholm_water_levels[[year - 1996]] <- stockholm_water_level
global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
global_avg_sea_level <- mean(mean_water_level_all_stations, na.rm = TRUE)
global_avg_sea_levels[year - 1996] <- global_avg_sea_level
} else {
cat("File for year", year, "not found. Skipping...\n")
}
}## File for year 2001 not found. Skipping...
## [1] -9.452900 3.475170 9.975784 9.191684 0.000000 2.850987 7.044084
## [8] 9.830954 8.804978 11.528670 19.242645 22.658903 19.638729 27.157361
## [15] 36.809982 33.308854 33.948932 33.369887 40.109891
There might be some repetition in the code and we could probably make it way easier but this is how it worked out in the end, but feel free to play around with it. Also we will iteratively plot the time series of the region to the global average.
years <- 1997:2015
yearly_stockholm_water_levels <- vector("numeric", length = 19)
for (year in 1997:2015) {
file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
if (file.exists(file_name)) {
dat <- nc_open(file_name)
water_level <- ncvar_get(dat, "waterlevel")
lat <- ncvar_get(dat, "lat")
lon <- ncvar_get(dat, "lon")
distances <- sqrt((lat - stockholm_lat)^2 + (lon - stockholm_lon)^2)
closest_station_idx <- which.min(distances)
stockholm_water_level <- water_level[closest_station_idx , ]
yearly_stockholm_water_levels[year - 1996] <- mean(stockholm_water_level, na.rm = TRUE)
global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
global_avg_sea_levels[year - 1996] <- mean(mean_water_level_all_stations, na.rm = TRUE)
} else {
cat("File for year", year, "not found. Skipping...\n")
}
}## File for year 2001 not found. Skipping...
manila_lat <- 14.5995
manila_lon <- 120.9842
yearly_manila_water_levels <- vector("numeric", length = 19)
for (year in 1997:2015) {
file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
if (file.exists(file_name)) {
dat <- nc_open(file_name)
water_level <- ncvar_get(dat, "waterlevel")
lat <- ncvar_get(dat, "lat")
lon <- ncvar_get(dat, "lon")
distances <- sqrt((lat - manila_lat)^2 + (lon - manila_lon)^2)
closest_station_idx <- which.min(distances)
manila_water_level <- water_level[closest_station_idx , ]
yearly_manila_water_levels[year - 1996] <- mean(manila_water_level, na.rm = TRUE)
global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
global_avg_sea_levels[year - 1996] <- mean(mean_water_level_all_stations, na.rm = TRUE)
} else {
cat("File for year", year, "not found. Skipping...\n")
}
}## File for year 2001 not found. Skipping...
######add manila to the plot
data <- data.frame(
Year = rep(years, 3),
Water_Level = c(global_avg_sea_levels, yearly_stockholm_water_levels, yearly_manila_water_levels),
Series = rep(c("Global Average", "Stockholm", "Manila"), each = length(years))
)
ggplot(data, aes(x = Year, y = Water_Level, color = Series)) +
geom_line(size = 1.2) +
labs(
title = "Global, Stockholm and Manila Water Level Changes (1997-2015)",
x = "Year",
y = "Water Level (mm)"
) +
scale_color_manual(values = c("blue", "red", "orange")) +
theme_minimal() +
theme(legend.title = element_blank())charlottetown_lat <- 46.2382
charlottetown_lon <- -63.1311
yearly_charlottetown_water_levels <- vector("numeric", length = 19)
for (year in 1997:2015) {
file_name <- paste0("hcc_obsclim_waterlevel_global_hourly_", year, ".nc")
if (file.exists(file_name)) {
dat <- nc_open(file_name)
water_level <- ncvar_get(dat, "waterlevel")
lat <- ncvar_get(dat, "lat")
lon <- ncvar_get(dat, "lon")
distances <- sqrt((lat - charlottetown_lat)^2 + (lon - charlottetown_lon)^2)
closest_station_idx <- which.min(distances)
charlottetown_water_level <- water_level[closest_station_idx , ]
yearly_charlottetown_water_levels[year - 1996] <- mean(charlottetown_water_level, na.rm = TRUE)
global_average_water_level <- rowMeans(water_level, na.rm = TRUE)
mean_water_level_all_stations <- apply(water_level, 2, mean, na.rm = TRUE)
global_avg_sea_levels[year - 1996] <- mean(mean_water_level_all_stations, na.rm = TRUE)
} else {
cat("File for year", year, "not found. Skipping...\n")
}
}## File for year 2001 not found. Skipping...
# plot
data <- data.frame(
Year = rep(years, 4),
Water_Level = c(global_avg_sea_levels, yearly_stockholm_water_levels, yearly_manila_water_levels, yearly_charlottetown_water_levels),
Series = rep(c("Global Average", "Stockholm", "Manila", "Charlottetown"), each = length(years))
)
ggplot(data, aes(x = Year, y = Water_Level, color = Series)) +
geom_line(size = 1.2) +
labs(
title = "Global, Stockholm, Manila, and Charlottetown Water Level Changes (1997-2015)",
x = "Year",
y = "Water Level (mm)"
) +
scale_color_manual(values = c("blue", "red", "orange", "green")) +
theme_minimal() +
theme(legend.title = element_blank())data <- data.frame(
year = rep(years, 4),
water_level = c(global_avg_sea_levels, yearly_stockholm_water_levels,
yearly_manila_water_levels, yearly_charlottetown_water_levels),
region = rep(c("Global Average", "Stockholm", "Manila", "Charlottetown"), each = length(years))
)
data$region <- factor(data$region, levels = c("Global Average", "Stockholm", "Manila", "Charlottetown"))
water_level_changes_plot.png <-ggplot(data, aes(x = year, y = water_level, color = region, linetype = region)) +
geom_line(size = 1.7) +
scale_color_manual(values = c("blue", "red", "orange", "green")) +
scale_linetype_manual(values = c("solid", "dashed", "dashed", "dashed")) +
labs(title = "Global and Regional Water Level Changes (1997-2015)",
subtitle = "Comparison of global average water level and regional locations",
x = "Year",
y = "Water Level (mm)") +
theme_minimal() +
theme(
legend.position = "top",
legend.title = element_blank(),
legend.text = element_text(size = 10),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
plot.title = element_text(size = 16, face = "bold"),
plot.subtitle = element_text(size = 12, face = "italic")
) +
scale_x_continuous(breaks = seq(1997, 2015, by = 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(min(data$water_level, na.rm = TRUE) - 10,
max(data$water_level, na.rm = TRUE) + 10))
print(water_level_changes_plot.png)Load the data and extract the water level information. Then calculate the mean of the waterlevels over time. This is done for the counterclimatic and observed values. By print() we can take a look at the data and gain some information on what is included, and most importantly we can see the names of the variables we would like to extract.
observed_hcc_relative <- nc_open("hcc_obsclim_waterlevel_global_monthly_1901_1978.nc")
print(observed_hcc_relative)## File hcc_obsclim_waterlevel_global_monthly_1901_1978.nc (NC_FORMAT_NETCDF4):
##
## 3 variables (excluding dimension variables):
## double lat[stations] (Contiguous storage)
## standard_name: latitude
## long_name: Latitude
## units: degrees_north
## axis: Y
## double lon[stations] (Contiguous storage)
## standard_name: longitude
## long_name: Longitude
## units: degrees_east
## axis: X
## double waterlevel[stations,time] (Contiguous storage)
## comment: This variable represents the monthly relative coastal water level as a combination of reconstructed geocentric water level data based on Dangendorf et al. (2019) and vertical land motion (VLM) based on Oelsmann et al. (2023). The data is aligned with the HCC hourly relative coastal water levels.
## coordinates: lon lat
## long_name: Relative Water Level
##
## 2 dimensions:
## stations Size:18719 (no dimvar)
## time Size:936
## long_name: Time
## axis: T
## units: days since 1900-01-16
## calendar: proleptic_gregorian
##
## 10 global attributes:
## Convention: CF-1.9
## featureType: timeSeries
## reference: https://doi.org/10.1038/s41558-019-0531-8
## references: The dataset is described in Treu et al. (2023) <https://doi.org/10.5194/essd-2023-112>. The source datasets are described in Dangendorf et al. (2019) <https://doi.org/10.1038/s41558-019-0531-8>, Muis et al. (2020) <https://doi.org/10.3389/fmars.2020.00263> and Oelsmann et al. (2023) <https://doi.org/10.5281/zenodo.8308347>
## institution: Potsdam Institute for Climate Impact Research (PIK)
## project: Inter-Sectoral Impact Model Intercomparison Project phase 3a (ISIMIP3a)
## contact: ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>
## isimip_id: 79f4a736-4f0d-4b24-a6af-88d99310b200
## grid_info: The irregular coastal grid points are distributed equidistantly along a smoothed global coastlines with a distance of 50 km globally and 10km in Europe with additional points added for tide gauge stations (described in Muis et all. 2019, https://doi.org/10.3389/fmars.2020.00263).
## title: Monthly coastal relative water levels
counterclim_hcc_rel <- nc_open("hcc_counterclim_waterlevel_global_monthly_1901_1978.nc")
print(counterclim_hcc_rel)## File hcc_counterclim_waterlevel_global_monthly_1901_1978.nc (NC_FORMAT_NETCDF4):
##
## 3 variables (excluding dimension variables):
## double lat[stations] (Contiguous storage)
## standard_name: latitude
## long_name: Latitude
## units: degrees_north
## axis: Y
## double lon[stations] (Contiguous storage)
## standard_name: longitude
## long_name: Longitude
## units: degrees_east
## axis: X
## double waterlevel[stations,time] (Contiguous storage)
## comment: This variable represents the counterfactual monthly relative coastal water level.It is derived by removing the quadratic trend since 1900 from the factual monthly relative coastal water level. The factual data is a combination of reconstructed geocentric water level data based on Dangendorf et al. (2019) and vertical land motion (VLM) based on Oelsmann et al. (2023).
## coordinates: lon lat
## long_name: Relative Water Level
##
## 2 dimensions:
## stations Size:18719 (no dimvar)
## time Size:936
## long_name: Time
## axis: T
## units: days since 1900-01-16
## calendar: proleptic_gregorian
##
## 10 global attributes:
## Convention: CF-1.9
## featureType: timeSeries
## reference: https://doi.org/10.1038/s41558-019-0531-8
## references: The dataset is described in Treu et al. (2023) <https://doi.org/10.5194/essd-2023-112>. The source datasets are described in Dangendorf et al. (2019) <https://doi.org/10.1038/s41558-019-0531-8>, Muis et al. (2020) <https://doi.org/10.3389/fmars.2020.00263> and Oelsmann et al. (2023) <https://doi.org/10.5281/zenodo.8308347>
## institution: Potsdam Institute for Climate Impact Research (PIK)
## project: Inter-Sectoral Impact Model Intercomparison Project phase 3a (ISIMIP3a)
## contact: ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>
## isimip_id: 65151046-bbdd-436b-b4b1-2809e39eb9a1
## grid_info: The irregular coastal grid points are distributed equidistantly along a smoothed global coastlines with a distance of 50 km globally and 10km in Europe with additional points added for tide gauge stations (described in Muis et all. 2019, https://doi.org/10.3389/fmars.2020.00263).
## title: Counterfactual monthly coastal relative water levels
Now we can take a first glimpse on the data we just extracted. The mean values can now for example be plotted on a world map to get a first idea of what the coastal water levels look like.
lat <- ncvar_get(counterclim_hcc_rel, "lat")
lon <- ncvar_get(counterclim_hcc_rel, "lon")
df <- data.frame(
lon = lon,
lat = lat,
mean_rel_obs = mean_rel_obs
)
world <- ne_countries(scale = "medium", returnclass = "sf")
ggplot() +
geom_sf(data = world, fill = "gray90") +
geom_point(data = df, aes(x = lon, y = lat, color = mean_rel_obs), size = 1) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, name = "Water Level (mm)") +
coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) +
labs(title = "Observed Coastal Water Levels (1901-1978 Mean)",
x = "Longitude", y = "Latitude") +
theme_minimal() +
theme(legend.position = "right")df2 <- data.frame(
lon = lon,
lat = lat,
mean_rel_counter = mean_rel_counter
)
ggplot() +
geom_sf(data = world, fill = "gray90") +
geom_point(data = df2, aes(x = lon, y = lat, color = mean_rel_counter), size = 1) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, name = "Water Level (mm)") +
coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) +
labs(title = "Counterclimatic Coastal Water Levels (1901-1978 Mean)",
x = "Longitude", y = "Latitude") +
theme_minimal() +
theme(legend.position = "right")In the next steps, we look at three different locations (Stockholm, Charlottetown & Manila) and plot the different sea level trends that can be seen in these regions.
First, let´s create a map with the locations highlighted:
stockholm_lat <- 59.3293
stockholm_lon <- 18.0686
stockholm_location <- data.frame(
lon = stockholm_lon,
lat = stockholm_lat,
name = "Stockholm"
)
charlottetown_lat <- 46.2382
charlottetown_lon <- -63.1311
charlottetown_location <- data.frame(
lon = charlottetown_lon,
lat = charlottetown_lat,
name = "Charlottetown"
)
manila_lat <- 14.5995
manila_lon <- 120.9842
manila_location <- data.frame(
lon = manila_lon,
lat = manila_lat,
name = "Manila"
)
world_plot <- ggplot() +
geom_sf(data = world, fill = "gray90", color = "white") +
geom_point(
data = charlottetown_location,
aes(x = lon, y = lat),
color = "red",
size = 2
) +
geom_text(
data = charlottetown_location,
aes(x = lon, y = lat, label = name),
hjust = -0.2, vjust = 0, size = 4, color = "red"
) +
geom_point(
data = stockholm_location,
aes(x = lon, y = lat),
color = "blue",
size = 2
) +
geom_text(
data = stockholm_location,
aes(x = lon, y = lat, label = name),
hjust = -0.2, vjust = 0, size = 4, color = "blue"
) +
geom_point(
data = manila_location,
aes(x = lon, y = lat),
color = "purple",
size = 2
) +
geom_text(
data = manila_location,
aes(x = lon, y = lat, label = name),
hjust = -0.2, vjust = 0, size = 4, color = "purple"
) +
labs(
title = "Locations on the World Map",
x = "Longitude",
y = "Latitude"
) +
coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)
)
print(world_plot)Next, we can plot time series for the single locations. This is done for the counterclimatic and observed data. On top of that we can calculate sea level trends for each data set and visualize this.
lat_diff <- abs(lat - stockholm_lat)
lon_diff <- abs(lon - stockholm_lon)
closest_station_index <- which.min(lat_diff + lon_diff)
#counterclimatic
stockholm_waterlevel_counter <- rel_waterlevel_counter[closest_station_index, ]
station_water_levels <- rel_waterlevel_counter[1, ]
time_period <- 1:length(station_water_levels)
#calculate the trend model
trend_model_counter_st <- lm(stockholm_waterlevel_counter ~ time_period)
slope_counter_st <- coef(trend_model_counter_st)[2]
intercept_counter_st <- coef(trend_model_counter_st)[1]
stockholm_data_counter <- data.frame(
Time = time_period,
WaterLevel = stockholm_waterlevel_counter
)
ggplot(stockholm_data_counter, aes(x = Time, y = WaterLevel)) +
geom_line(color = "lightgrey", size = 0.5) +
geom_abline(slope = coef(trend_model_counter_st)[2], intercept = coef(trend_model_counter_st)[1], color = "red", size = 0.5) +
labs(title = "Sea Level Trend for Stockholm (counterfactual data)",
x = "Time",
y = "Relative Sea Level (mm)") +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)
)stockholm_waterlevel_observed <- rel_waterlevel_observed[closest_station_index, ]
trend_model_observ_st <- lm(stockholm_waterlevel_observed ~ time_period)
slope_observ_st <- coef(trend_model_observ_st)[2]
intercept_observe_st <- coef(trend_model_observ_st)[1]
stockholm_data_observed <- data.frame(
Time = time_period,
WaterLevel = stockholm_waterlevel_observed
)
stockholm_data_observed <- data.frame(
Time = time_period,
WaterLevel = stockholm_waterlevel_observed
)
ggplot(stockholm_data_observed, aes(x = Time, y = WaterLevel)) +
geom_line(color = "lightgrey", size = 0.5) +
geom_abline(slope = coef(trend_model_observ_st)[2], intercept = coef(trend_model_observ_st)[1], color = "red", size = 0.5) +
labs(title = "Sea Level Trend for Stockholm (observed data)",
x = "Time",
y = "Relative Sea Level (mm)") +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)
)Now we can combine these two plots:
stockholm_combined_data <- data.frame(
Time = time_period,
WaterLevel = c(stockholm_waterlevel_observed, stockholm_waterlevel_counter),
Dataset = rep(c("Observed", "Counterfactual"), each = length(time_period))
)
plot_st <- ggplot(stockholm_combined_data, aes(x = Time, y = WaterLevel, color = Dataset, linetype = Dataset)) +
geom_line(size = 0.5) +
scale_color_manual(values = c("Observed" = "blue", "Counterfactual" = "orange")) +
scale_linetype_manual(values = c("Observed" = "dashed", "Counterfactual" = "solid")) +
geom_abline(slope = coef(trend_model_counter_st)[2], intercept = coef(trend_model_counter_st)[1], color = "orange", size = 0.9) +
geom_abline(slope = coef(trend_model_observ_st)[2], intercept = coef(trend_model_observ_st)[1], color = "blue", size = 0.9) +
labs(
title = "Observed vs Counterfactual Sea Levels in Stockholm",
x = "Time",
y = "Relative Sea Level (mm)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
plot(plot_st)lat_diff_charlottetown <- abs(lat - charlottetown_lat)
lon_diff_charlottetown <- abs(lon - charlottetown_lon)
closest_station_index_charlottetown <- which.min(lat_diff_charlottetown + lon_diff_charlottetown)
charlottetown_waterlevel_observed <- rel_waterlevel_observed[closest_station_index_charlottetown, ]
charlottetown_waterlevel_counter <- rel_waterlevel_counter[closest_station_index_charlottetown, ]
charlottetown_combined_data <- data.frame(
Time = time_period,
WaterLevel = c(charlottetown_waterlevel_observed, charlottetown_waterlevel_counter),
Dataset = rep(c("Observed", "Counterfactual"), each = length(time_period))
)
#observed
trend_model_observ_ch <- lm(charlottetown_waterlevel_observed ~ time_period)
slope_observ_ch <- coef(trend_model_observ_ch)[2]
intercept_observe_ch <- coef(trend_model_observ_ch)[1]
#counterfactuals
trend_model_count_ch <- lm(charlottetown_waterlevel_counter ~ time_period)
slope_count_ch <- coef(trend_model_count_ch)[2]
intercept_count_ch <- coef(trend_model_count_ch)[1]
#plot
plot_ch <- ggplot(charlottetown_combined_data, aes(x = Time, y = WaterLevel, color = Dataset, linetype = Dataset)) +
geom_line(size = 0.5) +
scale_color_manual(values = c("Observed" = "blue", "Counterfactual" = "orange")) +
scale_linetype_manual(values = c("Observed" = "dashed", "Counterfactual" = "solid")) +
geom_abline(slope = coef(trend_model_count_ch)[2], intercept = coef(trend_model_count_ch)[1], color = "orange", size = 0.9) +
geom_abline(slope = coef(trend_model_observ_ch)[2], intercept = coef(trend_model_observ_ch)[1], color = "blue", size = 0.9) +
labs(
title = "Observed vs Counterfactual Sea Levels in Stockholm",
x = "Time",
y = "Relative Sea Level (mm)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
plot(plot_ch)lat_diff_manila <- abs(lat - manila_lat)
lon_diff_manila <- abs(lon - manila_lon)
closest_station_index_manila <- which.min(lat_diff_manila + lon_diff_manila)
manila_waterlevel_observed <- rel_waterlevel_observed[closest_station_index_manila, ]
station_water_levels_manila <- rel_waterlevel_observed[1, ]
time_period <- 1:length(station_water_levels)
#observed
trend_model_observ_manila <- lm(manila_waterlevel_observed ~ time_period)
slope_observ_manila <- coef(trend_model_observ_manila)[2]
intercept_observe_manila <- coef(trend_model_observ_manila)[1]
manila_data_observed <- data.frame(
Time = time_period,
WaterLevel = manila_waterlevel_observed
)
#counterclimatic
manila_waterlevel_counter <- rel_waterlevel_counter[closest_station_index_manila, ]
trend_model_counter_manila <- lm(manila_waterlevel_counter ~ time_period)
slope_counter_manila <- coef(trend_model_counter_manila)[2]
intercept_counter_manila <- coef(trend_model_counter_manila)[1]
manila_data_counter <- data.frame(
Time = time_period,
WaterLevel = manila_waterlevel_counter
)
#plot
manila_combined_data <- data.frame(
Time = time_period,
WaterLevel = c(manila_waterlevel_observed, manila_waterlevel_counter),
Dataset = rep(c("Observed", "Counterfactual"), each = length(time_period))
)
plot_ma <-ggplot(manila_combined_data, aes(x = Time, y = WaterLevel, color = Dataset, linetype = Dataset)) +
geom_line(size = 0.5) +
scale_color_manual(values = c("Observed" = "blue", "Counterfactual" = "orange")) +
scale_linetype_manual(values = c("Observed" = "dashed", "Counterfactual" = "solid")) +
geom_abline(slope = coef(trend_model_counter_manila)[2], intercept = coef(trend_model_counter_manila )[1], color = "orange", size = 0.9) +
geom_abline(slope = coef(trend_model_observ_manila)[2], intercept = coef(trend_model_observ_manila)[1], color = "blue", size = 0.9) +
labs(
title = "Observed vs Counterfactual Sea Levels in Manila",
x = "Time",
y = "Relative Sea Level (mm)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
plot(plot_ma)For this I will show and adjust all the previous plots again. This is up to you, in what way you feel most comfortable doing this. Like this everything is together and it gives a better overview (if there are problems with the visualization it helps to export the final plot and look at it as a .png or pdf file):
plot_ma <- ggplot(manila_combined_data, aes(x = Time, y = WaterLevel, color = Dataset, linetype = Dataset)) +
geom_line(size = 0.5) +
scale_color_manual(values = c("Observed" = "blue", "Counterfactual" = "orange")) +
scale_linetype_manual(values = c("Observed" = "dashed", "Counterfactual" = "solid")) +
geom_abline(slope = coef(trend_model_counter_manila)[2], intercept = coef(trend_model_counter_manila)[1], color = "orange", size = 0.9) +
geom_abline(slope = coef(trend_model_observ_manila)[2], intercept = coef(trend_model_observ_manila)[1], color = "blue", size = 0.9) +
labs(
title = "Manila",
x = "Time",
y = "Relative Sea Level (mm)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "plain", color = "purple"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "none"
)
world <- ne_countries(scale = "medium", returnclass = "sf")
world_plot <- ggplot() +
geom_sf(data = world, fill = "gray90", color = "white") +
geom_point(
data = charlottetown_location,
aes(x = lon, y = lat),
color = "red",
size = 2
) +
geom_text(
data = charlottetown_location,
aes(x = lon, y = lat, label = name),
hjust = -0.2, vjust = 0, size = 4, color = "red"
) +
geom_point(
data = stockholm_location,
aes(x = lon, y = lat),
color = "blue",
size = 2
) +
geom_text(
data = stockholm_location,
aes(x = lon, y = lat, label = name),
hjust = -0.2, vjust = 0, size = 4, color = "blue"
) +
geom_point(
data = manila_location,
aes(x = lon, y = lat),
color = "purple",
size = 2
) +
geom_text(
data = manila_location,
aes(x = lon, y = lat, label = name),
hjust = -0.2, vjust = 0, size = 4, color = "purple"
) +
labs(
title = "Locations on World Map",
x = "Longitude",
y = "Latitude"
) +
coord_sf(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "plain"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10)
)
plot_ch <- ggplot(charlottetown_combined_data, aes(x = Time, y = WaterLevel, color = Dataset, linetype = Dataset)) +
geom_line(size = 0.5) +
scale_color_manual(values = c("Observed" = "blue", "Counterfactual" = "orange")) +
scale_linetype_manual(values = c("Observed" = "dashed", "Counterfactual" = "solid")) +
geom_abline(slope = coef(trend_model_count_ch)[2], intercept = coef(trend_model_count_ch)[1], color = "orange", size = 0.9) +
geom_abline(slope = coef(trend_model_observ_ch)[2], intercept = coef(trend_model_observ_ch)[1], color = "blue", size = 0.9) +
labs(
title = "Charlottetown",
x = "Time",
y = "Relative Sea Level (mm)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "plain", color = "red"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
plot_st <- ggplot(stockholm_combined_data, aes(x = Time, y = WaterLevel, color = Dataset, linetype = Dataset)) +
geom_line(size = 0.5) +
scale_color_manual(values = c("Observed" = "blue", "Counterfactual" = "orange")) +
scale_linetype_manual(values = c("Observed" = "dashed", "Counterfactual" = "solid")) +
geom_abline(slope = coef(trend_model_counter_st)[2], intercept = coef(trend_model_counter_st)[1], color = "orange", size = 0.9) +
geom_abline(slope = coef(trend_model_observ_st)[2], intercept = coef(trend_model_observ_st)[1], color = "blue", size = 0.9) +
labs(
title = "Stockholm",
x = "Time",
y = ""
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "plain", color = "blue"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "none"
)library(patchwork)
combined_plot <- (
(world_plot | plot_ch) /
(plot_ma | plot_st) +
plot_annotation(
title = "Sea Level Trends for different locations",
subtitle = "Observed vs Counterfactual Sea Levels for Selected Cities",
theme = theme(
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 14, hjust = 0.5)
)
)
)
print(combined_plot)